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Abstract 

We apply the coupled cluster method (CCM) in order to study the ground-state properties of the 
(unfrustrated) square-lattice and (frustrated) triangular-lattice spin-half Heisenberg antiferromag- 
nets in the presence of external magnetic fields. Approximate methods are difficult to apply to the 
triangular-lattice antiferromagnet because of frustration, and so, for example, the quantum Monte 
Carlo (QMC) method suffers from the "sign problem." Results for this model in the presence of 
magnetic field are rarer than those for the square-lattice system. Here we determine and solve 
the basic CCM equations by using the localised approximation scheme commonly referred to as 
the 'LSUBm' approximation scheme and we carry out high-order calculations by using intensive 
computational methods. We calculate the ground-state energy, the uniform susceptibility, the to- 
tal (lattice) magnetisation and the local (sublattice) magnetisations as a function of the magnetic 
field strength. Our results for the lattice magnetisation of the square-lattice case compare well 
to those results of QMC for all values of the applied external magnetic field. We find a value for 
magnetic susceptibility of x = 0.070 for the square-lattice antiferromagnet, which is also in agree- 
ment with the results of other approximate methods (e.g., x = 0.0669 via QMC). Our estimate 
for the range of the extent of the {M/Ms =)\ magnetisation plateau for the triangular-lattice 
antiferromagnet is 1.37 < A < 2.15, which is in good agreement with results of spin-wave theory 
(1.248 < A < 2.145) and exact diagonalisations (1.38 < A < 2.16). Our results therefore support 
those of exact diagonalisations that indicate that the plateau begins at a higher value of A than 
that suggested by spin-wave theory. The CCM value for the in-plane magnetic susceptibility per 
site is X = 0.065, which is below the result of the spin-wave theory (evaluated to order 1/5*) of 
XSWT = 0.0794. Higher order calculations are thus suggested for both SWT and CCM LSUBm 
calculations in order to determine the value of x for the triangular lattice conclusively. 
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I. INTRODUCTION 



Low-dimensional quantum magnets provide a difficult challenge to the theoretical physi- 
cist because of their strong quantum fluctuations and their complex dynamics [l, Q]. These 
effects lead to rich physics that include novel quantum phases, as well as quantum phase 
transitions between semi-classical magnetically ordered phases and magnetically disordered 
quantum phases, see, e.g., Ref. 

An interesting fleld of research is that of the behaviour of quantum magnetic systems in 
the presence of external magnetic flelds, see, e.g. Refs. 4, 5, 6, 7, ^. This topic has become 
more important by the discovery of exotic parts of the magnetisation curve of quantum 



antiferromagnets, such as plateaux and jumps 
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in the lattice magnetisation with respect to the externally applied fleld. Indeed, the presence 
of these plateaux and jumps may sometimes be linked purely to quantum effects because 
they are not observed in equivalent classical models at T = Q, Q, 21, 22|. Clearly, the 
behaviour of quantum magnetic materials in the presence of external magnetic flelds is an 
important aspect in their subsequent technological exploitation. Several methods such as 



quantum Monte Carlo method (QMC), fleld theories, exact diagonalisation of flnite systems, 
spin-wave techniques and strong-coupling approximation have been used 

mm m to 

study these systems. However, each method has its own speciflc limitations; for instance, 
the QMC is restricted (essentially) to unfrustrated systems because of the infamous 'sign 
problem.' 

In this article we focus on the behaviour of quantum antiferromagnets as they react to 
externally imposed magnetic flelds by a method of quantum many-body theory called the 



coupled cluster method (CCM) 
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32|. The CCM has been 



used previously in order to treat a wide range of strongly interacting quantum systems. 
In particular, the CCM is not restricted, in principle, by the spatial dimensionality of the 
problem or by the presence of competition between bonds, i.e., in frustrated quantum spin 
systems. A remarkable advance in the accuracy of the method for a localised approximation 



scheme called the LSUBm scheme has 
computer-algebraic implementations 



3een afforded by the use of "high-order" CCM via 
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Farnell and J Schulenburg 



29|. This computer code developed by DJJ 



is very flexible in terms of the range of underlying crystal- 



lographic lattice, spin quantum number, and types of Hamiltonian that may be studied. 
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Furthermore, recent advances to this code now allow "generalised expectation values" (with 
respect to one-spin and two-spin operators) and (separately) excited-state properties to be 
evaluated to high orders of approximation. Indeed, we employ the new code for the gener- 
alised expectation values to determine the lattice magnetisation and individual sublattice 
magnetisations of quantum antiferromagnets in external magnetic fields. 

The relevant Hamiltonian for an antiferromagnet in an external field is defined by 



H 



{id) 



(1) 



where the index i runs over all lattice sites on the lattice. The expression («, j) indicates a 
sum over all nearest-neighbour pairs, although each pair is counted once and once only. The 
strength of the applied external magnetic field is given by A. 

The quantum ground states at A = of all of the cases considered here are semi-classically 
ordered (albeit the classical order is reduced by quantum fluctuations) [2I. Classically, 
nearest-neighbours align in antiparallel directions for the bipartite antiferromagnets such 
as the antiferromagnet on the square lattice and at angles of 120° to each other for the 
Heisenberg antiferromagnet on the (tripartite) triangular lattice at A = 0. In the presence 
of an externally applied magnetic field (A > 0), the classical picture indicates that the spins 
will cant at various angles and that at a "saturation" value of A = Ag (square: A^ = 4; 
triangle; A^ = 4.5) all spins align with the field. The magnetisation saturates to a maximum 
value M = Mg at this point. 

However, we remark t 
net in a magnetic field ^ 



lat the behaviour of quantum spin-half square-lattice antiferromag- 

is (essentially) the same as that of the 
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classical model, a. 
spin- wave theory 



belt modified by quantum fluctuations. Second-order (and third-order) 
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37 



381 ] thus provides a good approximation to the behaviour of this 
model. Exact diagonalisations and QMC simulations [3, S] also provide good results for 
this case. Very recently, in Refs. 39|, |40|], the field dependence of the low-energy descriptors 



of this model (i.e., spin stiffness, spin- wave velocity, and magnetic susceptibility) have been 
investigated using exact diagonalisations and spin-wave theory. An excellent review of the 



properties of the spin-half square-lattice antiferromagnet is given by Ref. 



4l| 



By contrast, the behaviour o 



romagnet 



QflBfl 
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;he quantum case for spin-half triangular-lattice antifer- 
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43I ] is much different to that of the classical model. 
In particular, a magnetisation plateau is observed at M/Mg = ^ over a finite region of A. 



The range of this plateau has been estimated by spin-wave theory [lO|, to be given by 
1.248 < A < 2.145, whereas exact diagonahsations predict a region given by 

1.38 < A < 2.16. We note that the apphcation of the QMC method (leading to precise 
results for bipartite lattices) to the case of the triangular is severely limited by the "sign 
problem" due to frustration. The available spin-wave and exact-diagonalization data for the 
triangular lattice seem to be less accurate and complementary results are desirable. Fur- 
thermore, recent experimental evidence 1^ for the magnetic material Cs2CuBr4 suggests 
that a series of plateaux might exist at values of MjMg equal to 1/3, 1/2, 5/9 and 2/3. The 
authors of this article suggest that this might be due to unit cells of differing size for the 
different plateaux, e.g., each having an overall magnetisation of 1/2, and furthermore that 
theory has thus far only predicted the first of these at 1/3. However, the treatment of these 
possible higher plateau is beyond the scope of this article. 

The main goal of our paper is to explain how the CCM can be used to investigate the 
magnetisation process of quantum antiferromagnets and to provide detailed CCM results 
for the spin-half Heisenberg antiferromagnets on the square and the triangular lattices. The 
CCM has previously been applied with much success to the subject of quantum magnetic 
systems at zero temperature. The CCM provides accurate results even in the presence 
of very strong frustration. In particular, the use of computer-algebraic implementations 
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29( 1 of the CCM for quantum systems of infinite numbers of particles has been 



found to be very effective with respect to these spin-lattice problems. Here we present a 
brief description of the CCM formalism and its application via computational methods to 
the subject of quantum spin models. We then describe the application of the method to 
the spin-half Heisenberg model for the square and triangular lattices at zero temperature 
in the presence of an external magnetic field. We present our results and then discuss the 
conclusions of this research. 



II. THE COUPLED CLUSTER METHOD (CCM) 



As the CCM 
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29 



30, 



31 



las been discussed extensively elsewhere (see Refs. 
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24 



25 



26|, 



27 



32l|). we do not consider the methodology in depth here. In particular, the 



interested reader should note that the use of computer-algebraic implementations has been 



considered in Refs. 26, 



27 



28 



291 . However, it is still important to remark here that the 
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exact ket and bra ground-state energy eigenvectors, and (\^|, of a general many-body 
system described by a Hamiltonian if, are given by 

i7|*) = E,|*); {^\H = Eg{^\. (2) 

The ket and bra states are parametrised within the CCM as follows: 

|M/) = e^|$) ; S = J2SjCi, 

(§| = ($|^e-^ ; ~S = l + Y.SiCi. (3) 

One of the most important features of the CCM is that one uses a single model or reference 
state 1$) that is normalised. This, in turn, leads to a normalisation condition for the ground- 
state bra and ket wave functions ((\E'|\E') = ($|$) = 1). The model state is required to have 
the property of being a cyclic vector with respect to two well-defined Abelian subalgebras 
of multi-configurational creation operators {Cf} and their Hermitian-adjoint destruction 
counterparts {Cj = (C*/)^}- For spin systems the model state |$) typically can be chosen as 
an independent-spin product state and the corresponding operators {Cf} can be expressed 
as a product of a set of spin lowering operators, see below and for more details also Refs. 



26 



27 



28 



29|. 



The CCM formalism is exact in the limit of inclusion of all possible multi-spin cluster 
correlations within S and S*, although this is usually impossible to achieve practically. It 
is therefore necessary to utilise various approximation schemes within S and 5*. Here we 
use the localised LSUBm scheme, in which all multi-spin correlations over distinct locales 
on the lattice defined by m or fewer contiguous sites are retained. This approximation 
scheme has been successfully applied to determine the ground-state phases of quantum spin 
systems, see e.g. 



27 



29|. The CCM is a bi- variational formulation in which the bra and ket 
states are parametrised separately. This means that the ket and bra states are not explicitly 
constrained to be Hermitian conjugates. However, an advantage of this approach is that the 
Goldstone linked- cluster theorem is obeyed and so results may be found in the infinite-lattice 
limit N ^ oo from the outset. The important Helmann-Feyman theorem is also obeyed at 
all levels of approximation. The ket-state and bra-state equations are obtained using the 
following formulae, 

($|C7e-^iie^|$) =0, VJ ^ ; (4) 
($|5e-^[if,C+]e^|$) =0, V/t^O . (5) 
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FIG. 1: The model states used in the COM calculations for the Heisenberg model in an external 
magnetic field, a) The bipartite Lattices. Spins on the A and B sublattices make angles to the 
X-axis, b) The first model state for the triangular lattice (model state I). Spins on the A and B 
sublattices make angles a to the x-axis. Spins on the C sublattice point downwards, c) The second 
model state for the triangular lattice (model state II). Spins on the A and B sublattices point 
upwards. Spins on the C sublattice point downwards, d) The third model state for the triangular 
lattice (model state III). Spins on the A and B sublattices make angles a to the rr-axis. Spins on 
the C make an angle (3 to the x-axis. (Model state II is also a limiting case of model states I and 
III.) 



The methoc 
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in which Eos. @ and (Q are solved has been discussed extensively elsewhere 
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31 



32l |. However, we remark here that the computational 
method for solution of the CCM problem may be broken into three parts. The first task is, 
namely, to enumerate the fundamental set of CCM clusters for a given level of approxima- 
tion. Secondly, we must determine the ket-state equations in terms of the CCM ket-state 
correlation coefficients by pattern-matching those clusters Cy in the fundamental set to term 
in e~^He^ . Once we have determined the ket-state equations, the bra-state equations may 
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be determined directly. Finally, we solve the coupled CCM equations for the ket- and bra- 
state correlation coefficients, e.g., by using the Newton- Raphson method for the ket-state 
equations. Expectation values such as the lattice magnetisation may be obtained after we 
lave solved for both the ket and bra states. Again, we refer the interested reader to Refs. 



27 



28 



29| for more details of the practicalities of carrying out CCM calculations to 



high order. 

Here we use the classical ground states of these systems of the Heisenberg model in an 
external magnetic field as the model state. However, the magnitude of the characteristic 
canting angles in the quantum model (i.e., the angle between the local directions of the 
spins and the external magnetic field) may be different from the corresponding classical 
value. Hence, we do not choose the classical result for those angles. Indeed, we consider 
the angles as a free parameters in the CCM calculation, which has to be determined by 
minimisation of the CCM ground-state energy. 

The ground state of the classical system at zero external field (A = 0) has nearest neigh- 
bouring spins aligning in opposite directions for the bipartite lattices (e.g., the square lattice) 
and at angles of 120° to each other for the triangular lattice. Classically, the spins react to 
an external magnetic field by changing their alignment to that of the direction of the field. 
This is shown in Fig. [H For the bipartite lattices, the spins thus cant at an angle of 6 and 
vr — 6* to the x-axis, as is shown in Fig. [1^. By contrast, for the tripartite triangular lattice 
and related frustrated lattices one ought to distinguish between an applied field within the 
plane defined by the 120° planar state and a field perpendicular to this plane. Although on 



the classical level both cases are energetically equiva. 
fluctuations favour the planar configuration 
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ent 
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20 



21 



22l ] , thermal or quantum 



22| . Therefore in the present paper 



we restrict our considerations to planar states and a corresponding magnetic field applied 
within this plane. Following Ref. 15| we employ three different model states for the 
tripartite triangular lattice. The ffist such model state is one in which two spins on the A- 
and 5-sublattices point generally in the direction of the external magnetic field. However, 
they form angles a and vr — a to the X'clXlS, clS shown in the model state I of Fig. [T]d. The 
remaining spins on the C-sublattice point in a direction antiparallel to the applied external 
field. The second model state II of Fig. ^ for the triangular lattice has two spins on the A- 
and 5-sublattices that align completely with the external magnetic field and the remaining 
spins that align antiparallel to the external magnetic field. The final model state HI has two 



8 



FIG. 2: The bond directionality of the Heisenberg Hamiltonian after rotation of the local coordinate 
axes in the spin space. The directions of the bonds are indicated by the arrows placed on the square 
lattice. The two-site unit cell is also shown in dotted lines. 

spins on the A- and -B-sublattices that form an angle a to the x-axis and another spin on 
the C-sublattice that forms a (initially negative) angle of /? to the x-axis, as is also shown 
in Fig. [T]i. Model state II is clearly a limiting case of both model states, I and III. (For 
example, we obtain model state II from model state III by setting a = 7r/2 and f3 = — 7r/2.) 

In order to simplify the problem, we now rotate the local coordinate axes in the spin 
space so that all spins appear notationally to point in the downwards ^-direction. For a spin 
making an angle of 6 to the x-axis, the rotation of the local axes is given by, 

-s^sin(6') + s'=cos(6') 

-s^cos(^) - s"sin(^) . (6) 

The spins in the model state |$) now all appear to point downwards, i.e. |$) = |-- - iiii ■•■). 
The corresponding creation {Cf} are then given by {C/} = sf , sfs'j' , sfsj's'l , . . . , where 
the indices i,j,k,... denote arbitrary lattice sites. Furthermore, the Hamiltonian for the 
bipartite lattices in the rotated coordinate frame (i.e., with spins on the A sublattice making 
an angle 6 to the negative x-axis and spins on the B sublattice making an angle 6 to the 




positive shown in Fig. [1]) is now given by 



^ = E -7(1 + cos(2e))(44 + srsj) 



+ l(l-cos(2^))(.^T + ,rs;) 
-lsin(2^)(4.^^ + .-.^^)} 



+ Asin(^) Y.st- ^cosie) + srj 




+ 2^os(^)E(4, + ^rj- 



(7) 



We note that the sign in Eq. ([71) for those terms for s^s~^ and s^s~ for a bond going from i 
to j has an opposite sign for those same terms for a bond going from j to i. This is called 
a "bond directionality" and is indicated in the above equation by the arrow in the symbol 
{i j). An illustrative example of bond directionality in the Hamiltonian for the square- 
lattice case is shown in Fig. [2l We note also that runs over all A sublattice sites, is runs 
over all 5-sublattice sites, and i runs over all lattice sites. The translational symmetry of 
Eq. (171) compared to the original problem has also been reduced. We must include two sites 
in the unit cell, as is also shown in Fig. [2l 

Similar calculations may be carried out for the triangular lattice. We have three new 
Hamiltonians after rotation of the local spin axes of the spins for all three model states I, 
II, and III in Fig. [T](b-d) for the triangular lattice case such that all spins again appear to 
point downwards. The Hamiltonian for model state I, Fig. iH^b), for the triangular lattice is: 



H 



E 




+ 4(1 - cos(2a))(s+ s- + s- s+ ) 

- cos{2a)sl^sl^ + ^sin(2a)(s^^s+ + s^^s^ ) 




+ E 
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-cos(a)(s 



_|_ e- o-f 
^B,C ^C,A ' lB,C IC.A' 



- AEC + Asin(a)(EC+EC) 

ic iA ia 

- ^cos(a) + srj + ^cos(a) + srj , (8) 

M is 

where the sum {i^ is) goes from sublattice A to sublattice B (and with directionahty) . 
Note that {iB,c ic,A) indicates a sum that goes from sublattice B to sublattice C and 
sublattice C to sublattice A, respectively (and with directionality). A similar treatment 
may be carried out for the model state III, Fig. [I](d). Hence, if those spins on on the A and 
B sublattices make an angle a to the x-axis and those spins on the C sublattice make an 
angle /5 to the x-axis and employing the rotation of the local spin axes of Eq. ([6]), we find 
that, 



(ic^iA,B) 



1 

4 

+ cos(a 



^C ^A,B 

- ^sin(a - /3)(s^^s+ + s^^sr^ J 

+ E + ^^A<) + 

+ Asin(a)(E si + E < J + Asin(/3) E C 

+ ^cos(a){E« + 5rJ + E« + ^rJ} 

«A «S 

+ ^cos(/5)E« + .-), (9) 

^ ic 

where the sum {ic — > goes from sublattice C to sublattices A and i? (with directional- 
ity) and {lA, is) goes over each bond connecting the A and B sublattices, but counting each 
one once only (and without directionality). We note that we have three sites in the unit cell 
for all of the models states used for the triangular lattice antiferromagnet. 
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Note that in addition to the model states presented above, spin hquids such as valence- 
bond crystal states may be treated via the CCM is by using a dimerised or plaquette (etc.) 
as relevant model state. A corresponding matrix algebra 25| is then used with respect to 
this state. However, a simpler approach is now also available that relies on finding special 
solutions of the CCM equations for the Neel-type model states used here 32j. These allow 
us to treat via existing high-order formalism and computer code, for example, spontaneous 
symmetry breaking in the spin-half one-dimensional J1-J2 (Majumdar-Ghosh) model 32|. 
The CCM is thus not restricted purely to semi-classical systems. 

We consider the angles as free parameters in the CCM calculation. They are determined 
by direct minimisation of the CCM ground-state energy. This was achieved computationally 
at a given level of LSUBm approximation, and a minimum ground state energy with respect 
to these canting angles was also found computationally for a given fixed value of A. There was 
only one angle for the square-lattice antiferromagnet (and for model state I for the triangular 
lattice) and there were two such angles for model state III for the triangular lattice. The next 
value of A was then determined incrementally and the minimisation process of the energy 
with respect to the canting angles repeated. The fact that we had to minimise the ground- 
state energy with respect to such angles at each value of A made the CCM calculations 
much more costly in terms of computing time required than the equivalent situations at zero 
external magnetic field, which requires no such minimisation. Furthermore, we see that the 
Hamiltonians of Eqs. (1319]) do not conserve the quantity = J^i^i = 0, which is preserved 
for the square-lattice antiferromagnet at A = 0. For these reasons, CCM calculations in the 
presence of external magnetic fields are more challenging than their zero-field counterparts. 

A final point is that the inclusion of the CCM SUBl terms of form Si = Si-^^sf in 



the ground ket and bra states is also equivalent to a rotation of the local spin axes [23|. 
For example, for the spin-half system, we note that (s/')^|$)=0 and so we can prove that 
e'5'i|<|)) = 11^(1 +i5jjS^)|$). This produces a mixture of "up" and "down" spins at each site, 
which may be thought of (as may be seen from Eq. ([6]) above, for example) as the same as a 
rotation of local spin axes. Hence, we conclude that SUBl is equivalent to a rotation of the 



axes. Previous calculations for Heisenberg antiferromagnets in external magnetic fields [23 1 
made the explicit assumption that the correlation coefficients of the SUBl terms may be 
set to zero, and we make the same explicit assumption here. We minimise the ground-state 
energy explicitly with respect to the angles in our model state. Note that we go to much 
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higher orders of LSUBm approximation than those calculations presented in 23j. 

To investigate the magnetisation process in antiferromagnets we have to consider the 
total lattice magnetisation M along the direction of the magnetic field. This quantity (in 
the initial coordinate frame prior to rotation of the local spin axes) is defined by M = 
'^{Y.i^i) = Ei'^ll^') (s is the spin quantum number which is s = 1/2 throughout 

this paper). In the rotated coordinate frame (and in which all of the spins point appear 
"mathematically" to downwards), the lattice magnetisation for the bipartite lattices is now 
given by 

sin(^)^,~. cos(6'),^,~. , 

+ ^E(^ki + ^.^l^)^ (10) 

IB 

where, again, ia runs over all A sublattice sites, ib runs over all 5-sublattice sites, and i runs 
over all lattice sites. We are able to determine readily the lattice magnetisation once the 
ket- and bra-state equations have been solved for a given value of A. Furthermore, similar 
expressions to Eq. ( fTOl) may be obtained for the lattice magnetisation for the triangular 
lattice for model states 1, 11, 111, Fig. [Tl^b)-(d). We note that the magnetisation found on the 
three sublattices may become non-equivalent in a magnetic field for the triangular-lattice 
case. Indeed, for the triangular lattice, the expression for the lattice magnetisation aligned 
in the direction of the applied magnetic field on the individual sublattices (denoted. Ma, 
Mb, and Mc) in terms of the global axes prior to rotation of the local spin axes is given by 

MA,B,C = jr^^ E (^IC,.,cl^) ' (11) 

where the index ia runs over all A^^ sites on sublattice A, the index ib runs over all A^^ sites 
on sublattice 5, and the index iq runs over all Nq sites on sublattice C. Clearly, we see 
that N = Na + Nb + Nc and that M = {Ma + Mb + Mc) /3. 



III. RESULTS 



Now we present and discuss the results for the two models under consideration calculated 
by the CCM as illustrated above. We start with the spin-half square-lattice Heisenberg 
antiferromagnet. The ground-state energy in dependence of this model is shown in Fig. [31 
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0.5 1 1.5 2 2.5 3 3.5 4 



FIG. 3: Results for the ground-state energy per site Eg/N of the spin-half square-lattice Heisenberg 
antiferromagnet in dependence on an external magnetic field of strength A. Note that the curves 
for LSUB4, LSUB6, LSUB8 almost coincide. 
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FIG. 4: Results for the total lattice magnetisation M of the spin-half square-lattice Heisenberg 
antiferromagnet in the presence of an external magnetic field of strength A compared to results of 
QMC 0]. Note that the curves for LSUB4, LSUB6, LSUB8 almost coincide. 



14 



CT> 



0.5 



0.4 



0.3 



0.2 - 



0.1 




0.5 



1.5 



2 



2.5 



3.5 



FIG. 5: Results for the canting angle 9 /it obtained for the model state for the spin-half square- 
lattice Heisenberg antiferromagnet (see Fig. [1^) in the presence of an external magnetic field of 
strength A. Note that the curves for LSUB4, LSUB6, LSUB8 almost coincide. 



The CCM results converge rapid 
in previous CCM calculations 



y with increasing LSUBm level of approximation. As seen 
271], the ground-state energy in the limit of vanishing external 
field (A = 0) is approximated well. The interested reader is referred to Refs. 27|] for a more 
detailed discussion of these results. We also find that the exact result for the saturation field 
M = Ms at As = 4 is also reproduced. At this point the spins all lie in the direction of the 
external field. 

The results for the lattice magnetisation are shown in Fig. |H There is a considerable 
difference between the results for the spin-half quantum model and the classical straight- 
line behaviour (i.e., Mciassicai = |A). Clearly, this difference is because of quantum effects. 
It is also obvious from Fig. H] that the magnetisation of the quantum model is below that 
of the classical magnetisation in the region < A < A^. Again we note that the LSUBm 
results appear to converge with increasing m for all values of A. For example, the difference 
between the LSUB6 and LSUB8 results for the lattice magnetisation is less than 2 ■ 10"'^ for 
all values of A, and it is impossible to be detected by eye in Fig. HI From Fig. H] it is also 
evident that the CCM results for the lattice magnetisation are in excellent agreement with 



the results of QMC 



7|, which can be considered as the most accurate results available. 



In addition to the energy and the magnetisation we can also present results for the canting 
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angle 6 (cf. Fig. [T]) of the quantum model, see Fig. [51 Again, there is a noticeable difference 
between the values for the classical and the quantum angle. This difference first increases 
with A up to about A ~ 3.5. Beyond A ~ 3.5 the quantum angle very rapidly approaches 
the saturation value Qs = -n jl. 

In the next step the CCM results for the ground state energy and the lattice magnetisation 
in dependence on magnetic field can be used to calculate the uniform magnetic susceptibility, 
given by 

X.i^ = -1^. (12) 
Note that factor of \ in \^ is due to definition of M in the interval [0, 1]. Note further 



that we consider here x susceptibility per site 



4J] . For the concrete calculation of x we 



have used the second derivative of the energy. To check the accuracy for low fields we have 
also determined x numerically via direct determination from M by using We found 
that 2";^ and -j^-j^ agree to at least six decimal places of precision. 

The zero-field uniform susceptibility x(A — 0), the ground state energy, the sublattice 
magnetisation, the spin stiffness, and the spin-wave velocity constitute the fundamental 
parameter set that determines the low-energy physics of magnetic systems. The results 
for the ground state energy, the sublattice magnetisation, the spin stiffness for the square- 
lattice Heisenberg antiferromagnet at A = have been calculated by the CCM previously. 
The interested reader is referred to Refs. 27| for more details. However, CCM results 



for the susceptibility x were not determined by these earlier calculations. Here we find that 
;^=0.08596, 0.07915, 0.07650, 0.07498, and 0.07388 for the LSUB2, LSUB4, LSUB6, LSUB8, 
and LSUBIO approximations, respectively. Since the LSUBm approximation becomes exact 
for m oo, it is useful to extrapolate the "raw" LSUBm data to m ^ oo. Meanwhile 
there is much empirical experience how to extrapolate CCM LSUBm data for physical 



quantities such as the spin stiffness 



30, 



3l| and "generalised" susceptibilities 31| which 



are also related to a second derivative of the ground energy Eg. Hence, we use the same 
extrapolation rule for the zero-field uniform susceptibility that has previously been found 



to give good results for the spin stiffness and also for "generalised" susceptibilities [30, 
given by = cq + Ci/m + C2/m?- We see from Fig. [6] that this rule provides a good 

method of extrapolation of our data. The corresponding extrapolation then yields values 
for the susceptibility of x = 0.0700(6). (The number in brackets indicate the standard 
deviation.) This result is in reasonable agreement with data obtained by other methods, e.g 
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FIG. 6: LSUBm results for the zero-field uniform susceptibility x(A 0) for the spin-half square- 
lattice Heisenberg antiferromagnet (see Fig. [1^) with m = {2, 4, 6, 8, 10} and the polynomial fit 
according to x("^) = cq + ci/m + C2/m^. 



QMC (x = 0.0669(7)) |35|, series expansion (x = 0.0659(10)) [36^ 



(x = 0.05611) 3J], second-order spin- wave theory (x = 0.06426) : 
wave theory (x = 0.06291) |37|. 



inear spin-wave theory 
, and third-order spin- 



The field dependence of x is also of experimental interest, see e.g. [18|, |45|, |46|, |47|. We 
present LSUB4, LSUB6, and LSUB8 data for the field dependence of x in Fig- We 
note that the magnetisation divided by the applied external field is often considered in 
experimental studies. Hence, results for M/2A are given also in Fig. [71 For the sake of 
comparison, the classical value Xcias = 1/8 is also shown in this figure and we remark that 
this value is clearly independent of A. From Fig. [7| it is obvious that x and M/2X agree 
well with each other up to about A = 0.4 = A^/IO. The difference between results of the 
LSUB8 approximation and the classical result is about 4% at A = 0.4). However, these two 
sets of results begin to deviate significantly for larger A. Hence, the quantity M/2A is a 
good approximation for x for magnetic fields used in real experiments for systems with large 
saturation fields A^, and not for systems with low A^. We observe that x increases with A 
as we move away from the zero- field point, A = 0. Similar increases in x with the external 
field have been observed experimentally, e.g., for the quasi-two-dimensional antiferromagnet 



Ba2CuGe207 [45|. Moreover, these results are in agreement with recent results obtained 
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FIG. 7: Susceptibility x, see Eq. (fT^ . and the quotient M/2X in dependence on the magnetic field 
A the for the spin-half square-lattice Heisenberg antiferromagnet. Note that the M/2X curves for 
LSUB4, LSUB6, LSUB8 almost coincide. 



by exact diagonalisations, QMC simulations, and spin- wave theory 39|, |40|. As seen for 
these other methods, the susceptibihty is near the constant classical value for magnetic 
fields 1.5 ^ A ^ 3.5, although it starts rapidly to increase approaching the saturation field. 
Finally, weak oscillations seen for 1.5 ^ A ^ 3.5, although these are believed to be artefacts 
of CCM LSUBm approximation. We note that the number of oscillations increases are we 
increase the LSUBm approximation level, although their amplitude decreases markedly. In 
the limit, m — > cxd, it is expected that these oscillations will disappear entirely. 

We conclude from all of these results that the CCM provides precise results for the 
behaviour of the spin-half square-lattice quantum antiferromagnet in an external magnetic 
field. However, we see also from these results that the classical picture is essentially correct. 
Quantum mechanical effects modify, but do not change, the essential physics that occur in 
this unfrustrated quantum spin system. 

We now consider the spin-half antiferromagnet on the triangular lattice. However, the 
situation is more complicated here because we have three sublattices in this case. As dis- 
cussed above, we employ therefore the model states I, II, HI shown in Fig. [H^b-d). The 
computational effort of the CCM calculations presented here for the model state HI to very 
high orders is very great because we also need to find the minimum of the energy with re- 
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FIG. 8: Results for the ground-state energy per site Eg/N of the spin-half triangular-lattice 
Heisenberg antiferromagnet in the presence of an external magnetic field of strength A. Note that 
the results for LSUB4, LSUB6, LSUB8 are clearly converging rapidly for all values of A. 

spect to two canting angles, namely a and (3. The CCM calculation for the model state III 
in LSUB8 approximation was performed on a Beowulf cluster using 110 cores (Intel XEON 
3GHz CPU). On this computer the running time for one data point was approximately 2 
days. The CCM has been shown to be fully competitive with the results of other methods 
at the levels of approximation currently available to use using parallel computer methods 
(currently: a maximum of 1000 CPUs in parallel). The interested reader is referred, e.g., to 



Refs. 



26 



27 



28 



291] for detailed comparisons of CCM results to the best of other methods. 



The results for the ground-state energy are shown in Fig. [HI We note that the results for 
the model state with lowest energy are shown only as a function of A in Fig. [HI Thus, results 
of model state I only are presented for small values of the applied magnetic field strength 
A and results of model state III only are presented for higher values of A near to A^. The 
results of both model states coincide in the intermediate regime. Again, these LSUBm series 
of results are found to converge rapidly with increasingly levels of LSUBm approximation 
over all values of the external field parameter A. As may also be observed in Fig. [HI there 
is also a large reduction in the ground-state energy of the CCM results compared to the 
classical results for the energy (except in the trivial limit A — As = 4.5). 

The results for the total lattice magnetisation are shown in Fig. [9] The LSUBm results 
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are again seen to converge rapidly for increasing m. However, there is a radical departure 
from the classical straight-line behaviour (i.e. Mciassicai = |A) in this case. Thus, we find 
that the quantum model deviates from the linear relationship between M and A. The most 
prominent feature of our CCM results is the plateau in the M versus A curve at M/Mg = |. 
Note that the plateau corresponds to the "straight" part of the curve in the Eg{X) curve 
shown in Fig. [HI Note further that this plateau is well-known and has been found by other 
approximate methods 0, [?, ^ [9, Q, 1^. The ground state of the quantum system over the 
finite, non-zero range of A for the plateau region has ordering of the form shown in model 
state II of Fig. [H^c). Importantly, this is an example of when quantum fluctuations favour 
collinear ordering (so called 'order from disorder' phenomenon, see e.g. Refs. 48, 4^, 50|). 
This plateau state of model state II is observed only at a single point classically, namely, 
at A = 1.5. The classic al g round state is given by model state II in Fig. [T](c) only at this 



point, see also Ref. 
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22| . Indeed, states I, II and III are equivalent classically at the 



point A = 1.5. The values for the starting (Ai) and the end point (A2) of the plateau state 
calculated within different LSUBm approximations are shown in Table [B The most accurate 
values are provided by the LSUB8 approximation, namely, that Ai ~ 1.37 and A2 ~ 2.15. 
These results may therefore serve as the CCM estimate for the plateau width. We note that 
the results for Ai and A2 for even and odd values of m ought to converge to the same values 
in this limit. Our estimate for the range of the plateau is in also reasonable agreement with 
those results of spin- wave theory [l^ and exact diagonalisations j^j , which both predict a 
similar width for the plateau with respect to the applied external magnetic field. However, 
we note that spin-wave theory was carried out only to order 1/S for the triangular lattice 
antiferromagnet in an external field. We believe that higher orders than 1/S for spin-wave 
theory would provide better correspondence to those results of ED and CCM results cited 
here regarding the range of the plateau. The phenomenon of "order from disorder" in which 
quantum fluctuations tend to favour colinear states has studied extensively elsewhere, e.g.. 



Refs. 



48 



49| . We note that the plateau state (uud) is colinear in the present case, and so 



our results are another example of this phenomenon. We have shown here that quantum 
fluctuations stabilise the (uud) state over other states that classically would have had lower 
energy in the plateau region. 

We are able also to calculate the (sub)lattice magnetisation (i.e., with respect to the 
z-direction in the original unrotated spin axes) for the individual sub lattices, namely. Ma, 
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FIG. 9: Results for the total lattice magnetisation M of the spin-half triangular-lattice Heisenberg 
antiferromagnet in the presence of an external magnetic field of strength A. CCM results are 
compared to those results of exact diagonalisations [7|. The arrows illustrate the actual spin 
directions. We use model state I for A < Ai and we use model state III for A > A2 (see Fig. [TJ. 
Both model states give identical results within the plateau Ai < A < A2 • 



Mb and Mc given by Eq. ( jTTj) . by using the CCM and as a function of A. As far as we 
are aware, these quantities have never before been presented for this model. The results 
for Mai Mb and Mq are now presented in Fig. [TOl Once again, we see a radical shift in 
the quantum solution from the classical result. Interestingly, Mc appears to decrease before 
approaching the plateau at A = Ai, while Ma = Mb increase monotonically with A up to 
Ai. On the other hand, MaiMb decrease with magnetic field in the region A2 < A ^ 2.8 
above the plateau, while Mq increases monotonically with A up to A^. 

We discuss next the canting angles a and (3 in the model states I, II, III (see Fig. [TJ^a), 
(b), (c)) shown in Fig. [11] Note again that to the best of our knowledge data for the 
angles have not been presented previously by other authors. A strong difference between 
the results of the classical system and those results of the quantum system is again obvious, 
in particular, in the plateau region where in the quantum model a and (5 are constant 
but both angles change rapidly for the classical model. We see that the results for both 
a and (3 vary continuously, although not smoothly, for all values of A. There is no sudden 
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FIG. 10: Results for the sublattice magnetisation (7 = {A, i3,C}) on individual sublattices 
A and B (left) and C (right) of the spin-half triangular-lattice Heisenberg antiferromagnet in the 
presence of an external magnetic field of strength A. (Note that Ma = Mb for all A.) 
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FIG. 11: Results for the angle a/vr (left) and (3/tt (right) in the model state for the spin-half 
triangular-lattice Heisenberg antiferromagnet in the presence of an external magnetic field of 
strength A. The arrows illustrate the actual spin directions. 

discontinuity in the solution for the angles as was reported, e.g., for spiral phases of some 
frustrated quantum spin models. Note that above the plateau the angle a does not vary 
monotonously with field. Rather it first increases to a > 7r/2 reaching at maximum at about 
A ~ 3.2. Approaching the saturation then a rapidly decreases to a = 7r/2. As far as we 
aware, no such equivalent experimental results exist for the sublattice magnetizations or 
tilting angles. We recommend therefore that experimental investigations of these aspects of 
the magnetisation with external field also be carried out. 

For the zero zero-field uniform susceptibility ~^ 0), see Eq. (fT2|) . we obtain x=0.1139. 
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0.08568, 0.08200, and 0.07378 for the LSUBm approximation with m = 2,4,6, and 8. In 
addition, we can also calculate the individual response of the sublattices on the magnetic 
field, i.e. XA,Bfi = i^%r^- Due to the relation M = {Ma + Mb + Mc)/3 we have 
X = Xa+Xb+Xc- Again we can extrapolate the data for the susceptibilities to m — * oo using 
x{fn) = co+ci/m+C2/m^. The corresponding extrapolation then yields x = 0.065(23). (The 
number in brackets indicate the standard deviation.) We see from Fig.[T2]that this procedure 
is a reasonable method of extrapolation of the data for the triangular lattice, although it 
is not as good as for the square lattice. This is demonstrated by the magnitudes of the 
estimated standard deviations for the extrapolated values of x for the square and triangular 
lattices (of order approximately 10~^ and 10~^, respectively). We see from Fig. [T2l that the 
main contribution to x comes from the sublattices A and B. That is not surprising, since for 
the model state I, see Fig. [T]d, the direction of the magnetisation on the sublattice C is fixed, 
whereas the spins on sublattices A and B are rotated towards the field direction. Indeed, we 
find that Xa,b = 0.0245(54) and xc = 0.016(13) by extrapolating the susceptibilities on the 
different sublattices separately (see Fig. [T2l) . This analysis leads again to an overall value for 
x(= XA + XR + Xc) of Y = 0.065. We can compare this result with \_= 0.0794 obtained with 



spin-wave theory 



10 



42l ] . (We remark that this value of x in Ref ■ 42| was referred to as x± 



in this article and furthermore that it was defined per volume.) Although the magnitudes 
of X for the extrapolated CCM value and the spin-wave result agree, the difference between 
them is still obviously quite large. We believe that this difference might be attributed to a 
somewhat less reliable extrapolation (shown clearly in Fig. [T2|) than that presented for the 
square lattice above. However, we should note also that the spin-wave theory calculations 
of Ref. 42| were only ever carried out to order 1/S. (By contr ast, the spin- wave theory 



calculations for the square lattice were carried out to order 37|].) Hence, both higher 

order spin-wave results as well as higher order CCM-LSUBm results are recommended in 
order to establish a more accurate figure for x for the triangular-lattice case and, thus, to 
resolve this difference. 

Again we mention that the zero- field uniform susceptibility x(A 0), together with the 
ground state energy, the sublattice magnetisation, the spin stiffness, and the spin-wave 
velocity constitute the fundamental parameter set that determines the low-energy physics of 
magnetic systems. Corresponding CCM results for the ground state energy, the sublattice 
magnetisation, the spin stiffness for the triangular-lattice Heisenberg antiferromagnet at 
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FIG. 12: LSUBm results for the zero-field uniform susceptibilities xi^ ~^ 0) for the spin-half 
triangular-lattice Heisenberg antiferromagnet with m = {2, 4, 6, 8} and the polynomial fit according 
to x("i) = Co + ci/m -\- C2lrr? ■ 



A = can be found in Refs. [26|, |3C 

As for the square-lattice case above, we also present results at the LSUB4 and LSUB6 
levels of approximation for the field dependence of x Fig- [131 (Note that we have LSUB8 
data for x only for small fields due to the enormous computational effort of carrying out this 
calculation.) Again we compare x{^ with Mj (2A) which is often determined in experiments 
and also with the classical value Xcias = 1/9 that is independent of A. 

From Fig. [13] it is obvious that x ^-iid M /2\ agree well with each other up to about 
A = As/10 (the difference is about 7% at A = 0.45), but deviate significantly for larger A. 
As for the square lattice x grows with A starting from zero field up to the bottom of the 
plateau at Ai. In the plateau region x is zero indicating a finite excitation gap about the 
plateau ground state. Approaching the plateau from below or from above x{^) exhibits a 
sharp peak. Such peaks at the end of the plateau are indeed observed in experiments on an 
antiferromagnet on the triangular lattice, see e.g. Figs. 9 and 10 in Ref. [l^. Between the top 
of the plateau at A2 and the saturation at A^ we find a broad region where the susceptibility 
is small x ~ O-l- Approaching the saturation x again becomes large. The oscillations seen 
for A ~ 3.5 seem to be an artefact of CCM-LSUBm approximation. However, we expect 
again that the amplitude of oscillation will decrease with increasing approximation level and 
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FIG. 13: Susceptibility see Eq. (fT^ . and the quotient M/2X in dependence on the magnetic 
field A the for the spin-half triangular-lattice Heisenberg antiferromagnet. 

would disappear entirely in the limit m — oo. 



IV. CONCLUSIONS 

In this article we describe how the coupled cluster method (CCM) may be applied in 
order to calculate the behaviour of quantum antiferromagnetic systems in the presence of 
external magnetic fields. We have determined the ground-state energy, the total lattice 
magnetisation as well as sublattice magnetisations and the uniform susceptibility for the 
spin-half Heisenberg antiferromagnets on the square lattice and the triangular lattice by 
using the CCM to high orders of approximation. We showed that high-order CCM cal- 
culations give reasonable results for these quantities over all values of the magnetic field 
strength A for both lattices. For example, the CCM result for the lattice magnetisation for 
the square lattice compare well to QMC and spin-wave theory results for all values of the 
magnetic field strength. Our result for the uniform susceptibility of x = 0.070 for the square 
lattice is in reasonable agreement with those results of other methods (e.g., x = 0.0669(7) 
via QMC). Again, we believe that even closer agreement would occur with high orders of 
LSUBm approximation. 

CCM results presented here for the total lattice magnetisation for the triangular lattice 
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show the characteristic magnetisation plateau at M/Mg = | also seen in other studies 
0, 0, y, Q. The width of this plateau was estimated by us to be given by 1.37 ^ A < 2.15. 
This result was found to be in good agreement with results of spin-wave theory [l^ 
(1.248 < A < 2.145) and exact diagonalisations Q, 0, 0, 3 (1.38 < A < 2.16). Our results 
therefore support those of exact diagonalisations that indicate that the plateau begins at a 
higher value of A than that suggested by spin-wave theory. In addition, we provide results 
for sublattice magnetisations M^, Mb, and Mq evaluated on the individual sublattices A, 
B, and C of the triangular lattice that allows a better understanding of the magnetisation 
process of the triangular lattice. As far as we are aware, this is the first time that results for 
the individual sublattice magnetisations (and angles) have been presented. Our result for 
the longitudinal uniform low-field susceptibility x = 0.065 compares to the result of result 
of spin-wave theory [x = 0.0794), i.e. there is quite a large difference between the spin- wave 
and the COM result. Hence, higher order approximations for both SWT and COM LSUBm 
calculations and/or alternative approaches are recommended in order to obtain more 
reliable values for x for the triangular-lattice case. The susceptibility xW dependence 
on the magnetic field A shows for the triangular lattice characteristic sharp peaks at the 
bottom and the top of the plateau which may be used as indicators in experiments for a 
magnetisation plateau. 
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TABLE I: CCM results for the width of the magnetisation plateau for the spin-half Heisenberg 
antiferromagnet on the triangular lattice. 





Ai 


A2 


LSUB4 


1.312 


2.241 


LSUB5 


1.370 


2.030 


LSUB6 


1.357 


2.185 


LSUB7 


1.375 


2.105 


LSUB8 


1.370 


2.145 


SWT [lOj 


1.248 


2.145 


Exact Diagonalisations \7} 


1.38 


2.16 
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